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Abstract 

The recent observations of muon charge ratio up to about 10 TeV and of atmospheric neutrinos 
up to energies of about 400 TeV has triggered a renewed interest into the high energy interaction 
models and cosmic ray primary composition. A reviewed calculation of lepton spectra produced in 
cosmic ray induced extensive air showers is carried out with a primary cosmic ray spectrum that fits 
the latest direct measurements below the knee. In order to achieve this, we used a full Monte Carlo 
method to derive the inclusive differential spectra (yields) of muons, muon neutrinos and electron 
neutrinos at the surface for energies between 80 GeV and hundreds of PeV. Using these results the 
differential flux and the flavor ratios of leptons were calculated. The air shower simulator CORSIKA 
6.990 was used for showering and propagation of the secondary particles through the atmosphere, 
employing the established high energy hadronic interaction models SIBYLL 2.1, QGSJet-01 and 
QGSJet-II-03. We show that the performance of the interaction models allows makes it possible 
to predict the spectra within experimental uncertainties, while SIBYLL generally yields a higher 
flux at the surface than the QGSJet models. The calculation of the flavor and charge ratios has 
lead to inconsistent results, mainly influenced by the different representations of the K/tt ratio 
within the models. The influence of the knee of cosmic rays is reflected in the secondary spectra 
at energies between 100 and 200 TeV. Furthermore, we could quantify systematic uncertainties 
of atmospheric muon- and neutrino fluxes, associated to the models of the primary cosmic ray 
spectrum and the interaction models. For most recent parametrizations of the cosmic ray primary 
spectrum, atmospheric muons can be determined with an uncertainty smaller than ^3% of the 
average flux. Uncertainties of the muon- and electron neutrino fluxes can be calculated within an 
average error of and If 9%, respectively. 
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I. INTRODUCTION 



After a century from the discovery of cosmic rays, their origin is still the major quest 
in astrophysics. If ions are efficiently accelerated in diffusive shocks, supernova remnants 
in our Galaxy could be the major source of cosmic rays up to about 10 17 eV. In general, 
if hadronic particles are accelerated, a fraction of them must interact within their sources 
or in nearby molecular clouds to produce mesons, which eventually decay into high energy 
gamma rays and neutrinos with the energy spectrum ~ E~ 2 of the accelerated cosmic rays, 
while the rest propagates across the Galaxy until their detection on Earth. The observed all- 
particle cosmic ray spectrum can be generally described as a power law E' 1 , with 7 ~ 2.7 
up to ~ 3 x 10 15 eV (the so-called knee) and with 7 ~ 3.0 up to about 10 17 eV, above 
which energy, cosmic rays are believed to be of extragalactic origin. The interaction of these 
cosmic rays in the d ense Earth's atmosphere produce mesons and therefore, muons and 
neutrinos with a steep spectrum of ~ E~ 3 ' 7 pp. The search for fluxes of extra-terrestrial 
neutrinos rely on the fact that they should dominate at energy in excess of a few hundreds 
TeV over the atmospheric neutrino background. On the other hand, at this energy range 
the atmospheric neutrinos are mostly susceptible to large uncertainties due to cosmic ray 
spectral shape and composition above the knee, and to details of hadronic interaction models 
for the production of mesons in the atmosphere. A good understanding of the atmospheric 
neutrino flux is important for the identification of high energy cosmic neutrinos and therefore 
of the origin of the cosmic rays. 

Lepton fluxes in the atmosphere depend, to some extent, on cosmic ray composition 
(since this determines the proton/neutron, and therefore the charge/flavor ratio), and on 
the spectrum-weighted moments of the inclusive cross sections for the various hadronic 
processes, which depend on the high energy extension of hadronic interaction models. In 
particular, mesons produced in the extensive air showers (mostly pions and kaons), which 
can either decay into leptons or trigger further interactions with the atmosphere, shape the 
properties of secondary particles at the Earth's surface. The increasing relative contribution 
of kaons with energy is at the origin of the the observed increase of fi + / fi~ ratio in the 
TeV range [2H1] , because of the larger charge asymmetry from associated production of K + 
which is not compensated by a similar production of K~ in the forward regime [5]. The 
r ole of kaons is more important for neutrinos, since their energy is almost equally split 
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between the /x and the in the leptonic decay, while the pions' energy is mostly carried 
by the muons. Due to the steepness of the spectrum this effect is enhanced towards higher 
energies and kaons are the dominant source of neutrinos above a few hundreds GeV. The 
K/tt ratio also determines the response of muons and neutrinos to the seasonal variations 
in the stratospheric temperatures, which can be used to probe the contributions of heavy 
mesons in the extensive air showers [6HS] . Above a few 100 TeV the decay of heavy charmed 
mesons is expected to contribute to the lepton spectrum. Due to large uncertainties in charm 
production at large Feynman x, it is not known where such transition actually occurs, and 
it is model dependent. 



The atmospheric flux of + has been measured in a wide energy range from 1 
GeV to a few hundreds TeV. The IceCube Observatory reported the first determination 
of the high energy neutrino flux in the hundreds TeV region [9], [10] . In order to calculate 
the atmospheric neutrino intensities precisely, we need detailed information about (i) the 
primary cosmic-ray spectra at the top of the atmosphere, (ii) the hadronic interactions 
between cosmic rays and atmospheric nuclei, (iii) the propagation of cosmic ray particles 
inside the atmosphere, and (iv) the decay of the secondary particles. The comparisons 
between various calculations and with direct measurements makes it possible to assess how 
the contribution of experimental uncertainties in the primary spectrum and in the different 
hadronic models affects the atmospheric neutrino spectrum at the energy range that is 
relevant to the current neutrino telescopes. In this paper, the Monte Carlo calculation is 



done using the CORSIKA Extensive Air Shower Simulation code (http://www-ik.fzk.de/ 
corsika/). 



This paper is organized as follows: in Sec. [IT] we introduce the cosmic ray spectrum and 



composition, and in Sec. Ill the hadronic interaction models. In Sec. IV we describe the 
CORSIKA Monte Carlo calculation of lepton production in extensive air showers. In Sec. [V] 
we show the simulation results with comparisons with other calculations and direct exper- 
imental observations. First interaction models are benchmarked using muon observations, 
then the corresponding uncertainties on atmospheric neutrino spectra are determined. The 
effect on primary spectrum and composition is assessed as well. 



II. PRIMARY COMIC RAY SPECTRUM 



The observed primary cosmic ray flux covers a particle energy range from below 10 9 eV 
up to several 10 20 eV. In order to cover these 12 orders of magnitude in energy, a variety 
of different detection methods is used. Below ~ 100 TeV particle energy, air-borne and 
satellite experiments such as AMS [H], PAMELA [12], ATIC-2 p3], CREAM PUCE] and 
TRACER [16], directly measure particle energy and mass. Above ~ 100 TeV, the cosmic ray 
flux becomes too small and must be detected indirectly by large ground-based air-shower 
arrays (see Fig. [TJ. Different techniques such as the detection of secondary particles in 
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FIG. 1. All particle spectrum as measured by ground based arrays. The data are from |17H30| . 
The solid and dashed lines represent the power law models used as the parametrization of the 
primary cosmic ray flux for this work. Data compilation after |31j . 



an extensive air-shower or the measurement of air fluorescence are used. The fact that 
primary particles can only be observed indirectly by detecting portions of the induced air- 
shower makes it more difficult to measure the exact properties of the primary cosmic ray 
spectrum. A number of very large surface arrays have been built to observe cosmic rays up 
to about 10 20 eV. However, inconsistencies between the results of different experiments [28J, 



presumably originating in the deviation of the simulation from the data [32], show that a 
better understanding of hadronic interactions at ultra high energies is needed. 

The direct measurements provide the most unbiased results on cosmic ray spectrum and 
composition, and great progress was made the past years. In particular, recent observations 
by CREAM [HI 03] show an overall harder helium spectrum as compared to protons and, 
mainly, a flattening of their spectral index at about 230 GeV/nucleon. This result was 
confirmed by PAMELA [T2] and it is consistent with AMS results [llj. This flattening of 
the cosmic ray spectrum could be an effect of particle propagation through the Galaxy or 
of acceleration in their sources [33J, |3l] . 

The indirect measurements have provided the observation of cosmic ray spectrum at 
the knee region (which is at about 3 x 10 15 eV, where it steepens from ~ E~ 2 - 7 to ~ 
E~ 3 ) and up to the highest energies. Recent measurements by KASCADE-Grande [22] are 
a first observational hint that the knee in heavy particles shifts towards high energy, in 
general agreement with the notion of rigidity-dependence, although the result depends on 
the assumed hadronic interaction model. The highest energy cosmic rays above the so-called 
ankle at 3 x 10 18 eV are expected to be of extragalactic origin due to the high level of isotropy, 
which would not be present for Galactic sources [35J. 

In order to have a description of the expected atmospheric neutrino flux arising from 
cosmic ray interactions in the Earth's atmosphere, a careful parametrization of the cosmic ray 
composition is necessary. Some of the previous Monte-Carlo calculations reached neutrino 
energies up to 10 TeV e.g. Barr et al. [36] who uses a primary spectrum from Agrawal et al. 
[3T] or Honda et al. [3E] who use BESS [39J and AMS [llj cosmic ray data. Analytical 
calculations above the PeV region were performed by Sinegovsky et al. [10]. In particular 
the latter uses the model of the cosmic ray spectrum following Zatsepin and Sokolskaya [H] 
(ZS), who assumed three classes of Galactic sources. The first source class is the explosion 
of Supernovae into the interstellar medium, the second class is motivated by the explosion of 
supermassive stars into the local super-bubble and the third class explains the flux of nuclei 
below 300 GeV by Nova explosions. The ZS model provides a smooth transition from the 
all-particle spectrum measured in the direct experiments to that measured with extensive 
air showers, and it is compatible with the all-particle spectrum by KASCADE and GAMMA 
[12]. All considered models with a (rigidity-dependent) knee are motivated by the fact that 
both acceleration and propagation in models involving collisionless diffusion in magnetized 
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plasmas lead to the expectation of a rigidity-dependent cutoff for each individual component 
with a particle charge Z, E cut z oc Z |13H1S]- This can explain the steepening of the 
spectrum around the knee and can be taken into account in the modeling of the cosmic ray 
spectrum using a smoothed power law function as summarized in e.g. jHJ H7]. To effectively 
describe the all-particle spectrum of cosmic rays, five different primary mass groups, namely 
H, He, CNO, Mg-Si and Fe, are usually used to obtain a realistic representation, see e.g. 
[1HII19]. The individual spectra of the five components are summed up to get the all-particle 
spectrum. Recently, the PAMELA Collaboration has provided a new set of parameters for 
the proton and helium components of the first and third source class of the ZS model. 
These parameters are derived through a fit to their data p2]. The agreement to the data is 
significantly improved, thus in the following we use these updated parameters and refer to 
the model as (ZS/PAMELA). 

The poly-gonato model [SHU] describes the individual mass spectra up to the knee region 
fairly well, nevertheless the relatively steep dependence above the knee is not in agreement 
with the all-particle spectrum observations above about 10 17 eV. Primary particles in this 
energy range contribute to the production of leptons in extensive air showers in the 100 
TeV to PeV region. It is still disputed whether at this energy extragalactic cosmic rays 
can be considered as valid source class, or if a second Galactic component contributes to 
the primary spectrum between the knee and the ankle. In Hillas jSU] it is suggested that 
the primary cosmic ray spectrum is composed of three populations. The first population 
is associated to particles accelerated in Supernova Remnants with the knee indicating the 
cutoff. The second population (the so-called component B) accounts for the flux between 
the knee and the ankle and is associated with an unknown (or rather still debated) Galactic 
cosmic ray population. The third population is assumed to be of extragalactic origin and 
becomes dominant above about 10 18 eV. 

A recent ansatz including all three populations with up-to-date information on the com- 
position of the spectrum is presented in [3] . In this model each population is represented by 
5 mass groups with a simple power law and an exponential term representing the rigidity- 
cutoff. The spectrum for the first population is from the CREAM results extrapolated to 
a rigidity of about 4 PV to describe the knee, and the second population compensates for 
the all-particle disagreement between the knee and the ankle up to a rigidity of 30 PV. The 
extragalactic component is taken into account and two different composition scenarios are 
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investigated. The first one testing an all-proton approach (HGp) and the second one using 
a mixed composition (HGm) . 

In this paper, we will pursue the following approach: to obtain a realistic representation 
of the composition we follow the common approach to use the five-component model. Since 
the observed hardening of the primary spectra is incorporated into the HGp/HGm models, 
we prefer to use these as the default high energy model for this work. Both models are 
valid for energies above 10 TeV per particle. For a valid representation of the spectrum for 
energies below 10 TeV per nucleus, we extend the HG models to lower energies through a 
combination with the Gaisser-Honda 2002 (GH) [51] spectrum. The transitions between the 
two spectra are calculated taking the the high He contribution as: 1500 GeV for H, 2650 GeV 
for He, 7280 GeV for CNO, 45.8 TeV for MgAISi and 5050 GeV for Fe. In order to obtain a 
cross-over for the individual MgAISi fits of the two models, the parameters of the GH model 
are shifted within the given error boundaries, resulting in a different normalization constant 
KMg-si = 34.2 — 6 and a spectral index a>M g -si = 2.79 + 0.08. The air shower simulation is 
run for each component individually and it produces relative abundances of the secondary 
particles (muons and neutrinos). The results can then be weighted with different models of 
primary cosmic ray fluxes. 

The resulting all-particle spectra predicted by the various models are superimposed on 
air shower data in Fig. [TJ 

III. HADRONIC INTERACTION MODELS 

When a cosmic ray proton or nucleus enters the Earth's atmosphere, it can be treated 
as a projectile colliding with an air nucleus target. An interaction with their hadronic 
constituents results in the production of mesons and baryons which spread away from the 
original projectile track with transverse momentum p±. Unstable particles decay, while 
longer lived particles such as charged pions can undergo further interactions 7r =l= + Air — > 
anything, so representing a subsequent hadronic interaction with a different projectile-target 
configuration. The interplay of particle production, propagation and decay are the main 
ingredients of the atmospheric cascade. 

Air shower cascades are dominated by the soft component of the interaction, which rep- 
resents the hadronic cascades with small transverse momentum with respect to the shower 
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axis. Since no large momentum transfer is involved, the running coupling constant is too 
large for the application of ordinary perturbation theory [52] ■ Due to the lack of a self- 
contained theory for the soft phase space region, it is usually attempted to describe the 
physics with phenomenological approaches, such as the Regge theory. 

The air shower simulation code CORSIKA [53] includes various models for low- and high- 
energy interactions. For this study only the high energy interaction models are of interest, 
since the typical transition energy between high- and low-energy regime in CORSIKA occurs 
at the low energy boundary of this calculation at 80 GeV. We have selected the models 
SIBYLL 2.1 [M], QGSJet-Olc [55] with an additional heavy flavor (charmed) component 
and QGSJet-II-03 [56j, the successor of QGSJet-01. The restriction to these models is 
based on the acceptance in the air-shower community, the availability in CORSIKA and the 
computational time. While it would be interesting to test EPOS [57] it could not be done, 
since it demands approximately 60 times more calculation time than QGSJet-01. 

The Quark-Gluon-String models with minijet production (QGS Jet) are based on the phe- 
nomenological description of nuclear and hadronic collisions in Gribov's reggeon framework 
as multiple scattering processes [SSI ESI EH1 EH] • The individual scattering contributions, cor- 
responding to microscopic parton cascades, are described in terms of the exchange of "soft" 
and "semihard" pomerons. The Glauber approach is used for hadron- nucleus interactions. 
The original QGSJet model as well as QGSJet-II were especially designed for cosmic ray in- 
teractions, with the emphasis on the extrapolation to ultra-high energy cascades. QGSJet-II 
has been extended by the treatment of non-linear effects concerning very high energies and 
small impact parameters, where a large number of elementary scattering processes occurs 
and the underlying partonic cascades strongly overlap and interact with each other. An 
implication of this approach for hadron-nucleus and nucleus-nucleus interactions is that the 
non-linear screening effects are stronger in nuclear case, breaking the superposition picture. 
Additionally, the parameters of QGSJet-II are tuned according to more recent accelerator 
data corresponding to the state of 2006. Regarding the execution performance, QGSJet-II 
is 20 times slower than QGSJet-01. The heavy flavor generation of QGSJet-Olc results in 
the production of charmed hadrons, so the description of the prompt component of the at- 
mospheric lepton flux becomes possible. The model can handle charmed hadrons with the 
lowest mass only, i.e. neutral and charged D(D) mesons and A C (A C ) baryons [60J. These 
particles are explicitly transferred to the propagation code, given the chance for interaction 
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with air nuclei, which is rather uncertain due to the insufficiently known cross-sections. In 
CORSIKA, the decay is handled by Pythia 6.4 [HI] routines. Although, we have used 
this feature to explore the possibilities of the calculation of the prompt flux within a full 
Monte-Carlo approach, we emphasize that there are better and more reliable calculations 
available, exploring the prompt component [62TI66] . 

The underlying physical model in SIBYLL 2.1 [51] is the Dual Parton Model [52] for 
soft interactions with a mini-jet extension for the hard perturbative component [57]. Some 
features of the underlying code are borrowed from the Lund Algorithms contained in the 
Pythia 6 code. This model is explicitly optimized for air shower simulations, i.e. it im- 
plements extrapolation algorithms to ultra-high energies. Also, the program is efficiently 
designed to be called with a sequence of random collision energies and projectile-target con- 
figurations, in contrast to typical collider event generators, which expect to produce many 
events for the same configuration of incoming particles at the same center of mass energy. 

The mechanism of fragmentation/hadronization is similarly treated in SIBYLL and 
QGSJet. Both models describe the creation of new quarks and gluons in terms of one 
dimensional relativistic strings (color flux tubes), with one end attached to a valence (di-) 
quark from a projectile and a valence (di-) quark from the target, symbolizing the exchange 
of very soft gluons. The typical energy ("mass") density of such strings is in the order of 
k ~ 1 GeV/fm [61]. When the distance between the partons exceeds a certain critical value, 
the string breaks, creating a qq pair. The algorithm assigns each of these quarks to the open 
string ends, retaining the color confinement. Using this approach it is possible to explain 
the creation of new hadrons until the total energy stored in the original string has been 
dissipated as mass and forward momentum. 

Regarding the description of heavy nuclei collisions, both models use the Glauber ap- 
proach for the description of hadron-nucleus interactions. Major differences are related to 
the treatment of nucleus-nucleus interactions. SIBYLL implements the semi-superposition 
approach [HE], while QGSJet uses the quark-gluon-string picture. However, the calculation 
method of this work is based on averages of shower observables, so that the result will likely 
converge towards the super-position picture [lj 155]. 
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IV. OVERVIEW OF CALCULATION 



The central element of the calculation is the Extended- Air-Shower (EAS) simulator COR- 
SIKA version 6.990. We select a spherical detector to observe the theoretical flux without 
restricting the calculation to a certain type of detector. The simulation proceeds by injecting 
primary nuclei (H, He, C, Si and Fe) at the top of the atmosphere, such that propagation, 
hadronic interactions and particle decay are handled within the unmodified CORSIKA code. 
The temperature/density profile of the atmosphere is modeled by 5 exponential functions, 
with parameters fitted to the US Standard Atmosphere [69]. This static atmospheric model 
is widely used in neutrino flux calculations as a representation for the global atmosphere see 
[38] 1701 [71] , or [72] for an overview. Throughout the simulation, the Earth's curvature is 
taken into account as described in [73J. 

At lepton energies below hundreds of TeV the main sources of muons and neutrinos are 
decays of charged pions and various types of kaons (K ± , K s , K®). In the latter case two- and 
three-body decay modes are included. Using the hadronic generation counter in CORSIKA, 
it is only possible to identify leptons having a pion as mother particle or not, thus chained 
decays such as K — )■ ^ — )■ ji + are identified as pure pion decays. 

Due to constraints on the computational time, the simulation of the electro-magnetic 
component of EAS has not been activated, thus the channels having di-muon pairs in the 
final state via electromagnetic processes such as 7 — > fi + + \T are not included in this 
calculation. 

A. Calculation scheme and normalization 

When the secondary particles reach the sea level, they are binned according to their 
energy E, type p, mother particle m (K, n or charm), the high-energy interaction model 
Ai, the discrete energy of the primary nucleus E , the charge of the primary nucleus Z, the 
zenith angle of the primary nucleus 9 and the atmosphere A. This results in a database 
of yields y m ^ p (E p ,M,A,9,Z,E ). The yield is a differential inclusive energy spectrum 
(dN/dE) of a certain type of a secondary particle at the surface, as a result of an air shower 
cascade initiated by a number N of primary nuclei having several properties such as primary 
charge Z, zenith 9, interaction model A4 a.s.o. 
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The resulting flux of leptons of type p at the surface can be calculated as the discrete 
convolution 



%(e p , e,M,$ c ,A) = J2J2H W ( E ^ ®c{Zj)) x y mk -*p(E p , m, a, e, e „ z 3 ) , (i) 

m k Zj E ,i 

with the weight function 

w(E , ^c{Zj)) = s A ■ s N (E , QciZj)) . (2) 

The surface scaling factor sa compensates for the area normalization due to different refer- 
ence shells for flux of cosmic rays at the top of the atmosphere relative to the area of the 
virtual detector at the surface [72] : 

r E + h 

atm 1 m o /o\ 

s A = ~ 1.018 , (3) 

r E 

where r E is the Earth's radius and h atm is the atmosphere's height. The second factor 
sn{Eo^, Z, 9) compensates for the number of simulated showers N in the z-th primary energy 
bin AEq^, with respect to the physical flux in this energy bin according to some theoretical 
primary flux model $c : 

8 N (Eo, it 9 C {Z)) = f dE' $ C (E' , Z). (4) 

[Eo,i) JAE ,i 

The factor is calculated for each primary energy bin i and nucleus Z. In this work N(E = 
100 GeV) = 100, 000, extending with a power law behavior (~ E^ 1 ) up to N(E = 
100 EeV) = 40,000. This corresponds to a differential spectral index of 7 = 1.05. This 
approach allows us to efficiently re-weight the simulated dataset according to any primary 
model without the recomputation of the air shower database. 



B. Approximations 

For secondaries above 80 GeV the 3D-effect, the east-west effect and the up-down asym- 
metry have a negligible contribution [36j EH EI]- Also, for primaries above hundreds of 
GeV the influence of the geomagnetic cutoff and the Earth's magnetic field is below the 
simulation accuracy. Therefore, we make following approximations: 

1. the flux of primary cosmic rays at the top of the atmosphere is isotropic, 
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2. the error due to joint usage of the US Standard atmosphere is smaller than other 
systematic uncertainties, 

3. the effect of Earth propagation is neglected, thus the flux of neutrinos is assumed to 
up-down symmetric up to ~ PeV energies (see e.g. [7_H ES] for the estimation of this 
effect on the event rates of neutrino telescopes) 

4. the flux is azimuth-symmetric, 

5. neutrino oscillations have no significant effect, 

6. the lateral distribution of the EAS is not taken into account and particles are consid- 
ered to travel on a line trajectory, even if the full tree-dimensional cascade is simulated. 

Since in CORSIKA the shower has a fully three-dimensional shape, the last approximation 
is motivated by the disregard of the lateral distribution. The total number of simulated 
showers is 787,161,600. 

V. RESULTS 

The structure of this chapter is as follows: in the first part the particle charge and 
flavor ratios are calculated, reflecting some physics assumptions of the employed interaction 
models and features of the air shower code. These values are mainly sensitive to the pion to 
kaon ratio in the air shower. The differences between the interaction models are dominant 
compared to the variation due to the primary cosmic ray flux. Therefore, the calculation of 
the conventional ratios has been carried out, using the combined Gaisser- Honda and Hillas- 
Gaisser model with protons only in the extragalactic component (cHGp). If the zenith range 
is specified as vertical, then it is restricted to cos(6') < 0.25 and cos(#) > 0.75 for horizontal, 
respectively. Generally, neutrino fluxes are given as averages over all zenith angles and muon 
fluxes for the vertical direction only. 

In the second part, we present the calculated conventional fluxes of atmospheric muons, 
muon neutrinos and electron neutrinos. Here, we again use the cHGp model to parametrize 
the primary cosmic ray flux and composition. 

In the last part, we study the influence of the knee of cosmic rays by variation of the 
primary CR flux model for a given interaction model. Additionally, we show the prompt 
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prediction based on the implementation of QGSJet-Olc for different primary flux models. 
A. Interaction model performance 

The ratios of neutrino flavor and of muon charge represent the interaction model perfor- 
mance regarding the inclusive spectra of pions and kaons in the atmospheric cascade, where 
kinematical propagation effects, a realistic atmosphere and air composition are taken into 
account. There is also sensitivity to the primary composition, since a higher contribution 
of neutrons in the all-nucleon flux lowers the fraction of positively charged particles in the 
shower. The statistical uncertainty of the Monte-Carlo method limits the reasonable energy 
range to < 100 TeV. The errors are purely statistical. 

1. Lepton ratios at the surface 

Figure [2] shows ratios of neutrino flavor and of muon charge as they would be observed 
by a surface detector. In the upper left panel, the muon neutrino to antineutrino ratios are 
drawn together with some reference calculations. HKKM 2011 [38], which uses a modified 
version of DMPJET-III [78 J for the high energy part, performs close to the cascade equation 
approach by Sinegovsky et al. [TH] , employing the parametrization of nuclear interactions by 
Kimel-Mokhov and the ZS primary flux model. The Bartol 2004 [36J calculation reflects the 
properties of the interaction model TARGET [79], which explicitly includes the associated 
kaon production channel p + anything — > A + K + [37J |5U [79]. The similarities between 
SIBYLL and TARGET become apparent due to the higher kaon charge ratio and thus the 
increased over flux, due to the similarly treated associated kaon production. Both 
QGSJet models fall below all other expectations because the charge ratios (bottom panels 
of Fig. [2]) of pions and kaons in the air shower are constant over a wide energy range, showing 
that no associated kaon production is included. 

In analogy to the i^/F^-ratio, the muon charge ratio for vertical zenith directions is shown 
in the top right panel of Fig. |2j As expected from the muon neutrino ratios, the muon charge 
ratio from SIBYLL suffers from the overestimation of K + . However, the shape of the curve 
reproduces the features observed by experiments. It is therefore presumably a matter of 
scale between the inclusive pion and kaon spectra in the p N-/n iV-interaction or the ratio 
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FIG. 2. Top left and middle panels: All angle averaged ratios of neutrinos, compared to calculations 
[36] . [38] . [76] and [10]. Top right panel: Muon charge ratio for vertical muons, compared to data by 
MINOS [2 J and L3+C [77J. Bottom panels: Charge ratio of pions (left) and kaons (right) derived 
from surface muons (full markers) and muon neutrinos (hollow markers), respectively. 
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of the a-jr-Air to the (Tk-ait cross-sections. Due to the flat meson charge ratios and also the 
low kaon charge ratio, both QGSJet models do not describe the shape and normalization 
according to L3+C and MINOS data correctly. 

The middle left panel shows the i/ e /F e -ratio. Again, the difference between the interaction 
models which have an increased K + production can be observed due to the production 
channel p — > AK + and the three-body decay K + — > 7[°e ± v e {v e ). Within SIBYLL this 
process has a higher contribution (see bottom panels). 

The neutrino flavor ratio in the middle right panel is sensitive to the occurrence of the 
3-body kaon decay process K^ ev [81)] , which is the dominant source of conventional electron 
neutrinos for the observed energy. Our calculation is close to HKKM 2011 for all three 
interaction models. From [ID] the SIBYLL + GH calculation is used. The falling ratio for 
QGSJet-01 at E u > 10 TeV is from the additional v e flux, coming from decays of charmed 
hadrons. 

The bottom panels show the charge ratio of the mother mesons of muons and neutrinos at 
the surface. In CORSIKA, secondary neutrinos are treated as final state particles without 
further interactions. The spread of the muon charge originating from pions for energies 
below 1 TeV only occurs for horizontal zenith angles and thus it can be explained through 
the occurrence of chained decays of K + — > n + + anything — > fi + + and muon decay 
\i — > e + v e {pe) + ^tii^fi)- In the right panel, the kaon charge ratio is nearly constant for 
all energies and the QGSJet models (no associated production). In SIBYLL the enhancement 
due to associated kaon production is apparent as stated above. This behavior suggests that 
the A cross-section is too large or their spectrum too hard. This effect enhances the total 
kaon multiplicity in the cascade, leading to a higher neutrino flux at the surface. 

2. Mesonic origin of leptons 

In Fig. [3] the fractions of muons and neutrinos are shown with respect to their mother 
particle of arbitrary charge. As expected, SIBYLL has a higher contribution of kaons at 
all energies. For energies above E v > 100 TeV the source of conventional neutrinos are 
nearly exclusively kaons. Superimposed in the top and middle panels are curves, which are 
obtained by combining the approximative analytical formulas (6.4) and (7.5) from pQ with 
experimentally fitted parameters from MINOS [2] (labeled as MINOS pi-K in the figure) 

16 




0.7 
0.6 
0.5 



'□□□□□□□□□□□aOi ^ Jji 



ji-* fl, SIBYLL2.1 • K- n, SIBYLL2.1 

n - n, QGSJET-II ■ K-> ft QGSJET-II 

q ^ — Jt -* |x, Minos pi-K - -K -* (x, Minos pi-K 

0.3 



0.2 

0.l|»v.*"-'""rt 





.v*'v«.-***++±T TI 



; 



horizontal 



10 2 



10 3 



10 4 



10 5 10 6 
E M / GeV 



e 1 0.9 



0.8 r 




'□□□□□ nDnp n n $ tfr 



0.7 r ^"oOooS^-oOqUooO*^^^! 

0.6 r O It-* fi, SIBYLL 2.1 •«-»[*, SIBYLL 2.1 

5 - □ it - |t, QGSJET-II ■ K-*n, QGSJET-II 
p — Jt-*m Minos pi-K - -K -* |i, Minos pi-K 

04 . ... r? >-v^,: s .,!! M . 

02 - A TI i 

•••,,-'.«■ t 

01 iV"'"'"" vertical 



10 2 



10 3 



10 4 



10 5 10 6 

/ GeV 



0.5- 



© 



-*MA t 



□ 7t->n. QGSJET01C 

0.5 k -> 11, qgsjetoic 

* charm -» n, QGSJET01C 



if 

© 



8« 



* charm -> \i, QGSJET01 c ik Slri t 

■■■••Vil+f yi, 

vertical * * "* t+t Ik , ■ 

10 2 10 3 10 4 10 5 10 6 10 7 10 s 

E^, / GeV 



i 1 

o 

J? 0.9 
0.8h 

V ^ 

0.7: 



o.6r°o D a 

PI • 



o it - v p , SIBYLL 2.1 i K-* v |( , SIBYLL 2.U 
□ it - v„, QGSJET-II ■ K - v„, QGSJET-llj 




10° 
EJ GeV 



• .■■ ■ i 



o it -v fl SIBYLL 2.1 . K-»v |( , SIBYLL 2.1 i 
□ it - v , QGSJET-II ■ K - v , QGSJET-II 4 
— it -* v , analytic _-K^v analyic " 




^0.5 
© 



o jt -» v^, QGSJET01 c 

• K — » v^, QGSJETOIc 

* charm -*y QGSJETOIc J 



vertical 



10 2 10 3 10 4 10 5 10 6 10 7 



10° 

E v / GeV 



FIG. 3. Fractions of pions and kaons to the total flux for vertical and horizontal muons (left hand 
side) and neutrinos (right hand side), respectively. The solid curves are obtained through the 
approximative formulas from [1], a detailed explanation is in the text. The four bottom panels 
separately show the corresponding figures obtained with QGSJetOlc including charmed hadrons. 
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or the SIBYLL-2.1 at 10 TeV z-factors from [81J (labeled as analytic in the figure). In the 
muon channel (left panels) the MINOS parametrization suggests that the fractional meson 
contributions of the air shower are better described by QGSJet at energies below 300 GeV, 
but at higher energies they approach the SIBYLL model. 

The four bottom panels show QGSJetOlc predictions including the prompt component. 
For the vertical incident direction the conventional flux is suppressed by the steeper atmo- 
spheric pressure gradient, such that longer lived mesons prefer the interaction with air nuclei 
prior the decay. The very short lived .D's and A c 's promptly decay into leptons, which carry 
a large fraction of the original energy of the nucleus. The model in QGSJetOlc predicts that 
the purely prompt flux can be observed above 10 PeV for muons or neutrinos irrespective 
the direction of incidence. 



B. Differential lepton fluxes at the surface 

In contrast to the previously presented ratios, the total muon neutrino intensity can 
directly be observed by recent neutrino detectors, such as the IceCube Observatory and 
ANTARES. However, the uncertainty of the measurements is to some extent dependent on 
the assumed model of hadronic interactions, the atmosphere and the primary flux during 
the data analysis. Beside the cHGp primary flux parametrization usage throughout the 
calculations of the next subsections, we employ the GH spectrum when it is possible to 



compare the result of this study with measurements. In Sec. |V C| the influence of the 
primary model on the fluxes is studied explicitly. The neutrino fluxes are averaged over the 
azimuth and zenith angles, while the muon fluxes are presented for the vertical direction 
only. The prompt component is turned on in all simulations with QGSJetOlc. 

The differential flux of conventional + is presented in Fig. [4j compared to experi- 
mental data from AMANDA II and IceCube with 40-string configuration. 

The flux obtained with the two QGSJet models is is comparable with Amanda and Ice- 
Cube data in the experimentally observed range. As discussed in the previous section, 
SIBYLL has a higher fraction of kaons leading to a higher muon neutrino flux. In connec- 
tion with the harder primary spectrum compared to the previous calculations, it seems to 
overestimate the flux with respect to the IceCube data. In the region of hundreds of GeV, 
where the GH primary model is valid, the Monte-Carlo calculation agrees well with the 
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FIG. 4. Flux of conventional Vn + averaged over zenith and azimuth angles, compared to 
experimental data from Icecube in 40-string instrumentation [9j and Amanda II [10] . The bottom 
pane shows the full energy range for muon and electron neutrinos, calculated using the cHGp model 
and compared to calculations by |40j , [38] and [36] , providing the visibility of the influence of knee 
of CR in muon neutrino spectra. 



semi-analytic calculations from Sinegovsky et al. [40J. 



19 



The calculation of the flux of electron neutrinos is depicted in the bottom pane of in Fig. 
[4] for the full energy range of the simulation. Due to the lack of experimental observations 
the results can only be compared to other calculations. The flux calculated with SIBYLL 
or either version of QGSJet agrees with other calculations. The significant contribution of 
prompt neutrinos practically eliminates the effect of the knee on the spectral shape with 
QGSJetOlc. Therefore, every single electron neutrino detected with E v > 1 PeV has a very 
high probability to have a prompt or astrophysical origin. 

In Fig. [5] the muon flux is compared to data and calculations. The spread between the 
different interaction models decreases due to the minor contribution of kaons. In the range 
of the pion's critical energy (e n ~ 115 GeV), the BESS and L3+C data strongly constraints 
the validity of the calculation to either combination of SIBYLL-2.1 or QGSJet-II using the 
cHGp spectrum. The QGSJet-Olc simulation yields a too low muon intensity with either 
primary spectra below TeV energies. 

C. Variation of the primary cosmic ray flux 

Figure [6] shows the muon charge ratio, calculated with SIBYLL-2.1 using different primary 
flux models. Although SIBYLL's prediction is too high, it does well reproduce the transition 
between the pion charge dominated (E < e n ) and the kaon charge dominated (E > ex) 
regions of the charge ratio. As pointed out in [5], the variation of the primary flux model 
leads to the variation of the slope of this transition, which is steepest for the ZS spectrum 
and flattest for the two HG models. 

In Fig. [7j we have calculated the surface fluxes assuming different primary spectra and 
compositions for the primary cosmic ray flux with respect to a baseline spectrum. To 
emphasize the differences of this calculation in connection with previously published primary 
cosmic ray flux models, we have selected GH (2002) as the baseline. The results are similar 
for QGSJet-Olc. The shape of these curves does not change when using SIBYLL-2.1, but 
the features are shifted roughly a factor ~ 2 towards higher energies in the case of muon 
neutrinos and a factor of ~ 4 in case of muons, i.e. the ratio $^(A^ = cHGp)/<l ) M (A / l = GH) 
crosses unity at 800 TeV instead of 200 TeV. 

The poly-gonato model yields the lowest flux, falling below all other models above 500 
GeV. This model is designed with the goal to describe the cosmic ray flux below the knee 
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FIG. 5. Surface flux of vertical (on the left pane for cos(#) = 1.0) fi + + compared to measure- 
ments and calculations. The cascade equation calculations are from [81J. The experimental data 
are from [39" t 177 1 I82H86] and datapoints for [82H85] have been taken from [5T] . 

and at the knee. Above, the spectrum is to steep and does not agree with data (see Fig. [T]). 
It is therefore not suited to accurately describe effect of the knee on atmospheric leptons. 

The Zatsepin-Sokolskaya (PAMELA parameters) model agrees with several indirect mea- 
surements at energies close to the knee and with direct PAMELA measurements in the 
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FIG. 6. The charge ratio, deduced from the vertical muon flux, using SIBYLL-2.1 and different 
primary models. The MINOS 7r-K model and the data are from [2] and for L3+C [77] respectively. 
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proton and helium component. Using this model the lepton fluxes show a significant kink 
at tens of TeV, originating from the transition of the first (SN) to the second (SN into 
super-bubble) source class. This transition leads to a variation of the lepton fluxes in the 
order of 20 - 30%. 

The two Hillas-Gaisser models (cHGp and cHGm) incorporate the hardest spectrum, and 
thus lead to the highest fluxes at lepton energies above several TeV. The hypothetic second 
Galactic component plays an important role at the knee, being the source of atmospheric 
leptons at knee energies. Due to the overall good agreement of these models with the all- 
particle primary spectrum we favor the version with protons in the extragalactic component 
for all our neutrino calculations. 

The effect of the knee of cosmic rays in the primary spectrum is reflected in a similar 
shape for muons and muon neutrinos. However, the logarithmic abscissa does not represent 
well the differences in energy. The spectral index begins to change at several TeV and falls 
below the knee-less hypothesis (GH spectrum) in the range 100 - 200 TeV in the case of the 
ZS and the cHGp/m models. 

D. Charm in QGSJet-Olc 

We have used the implementation of charm hadron production in QGSJet-Olc to assess 
the influence of the primary spectrum and composition on the prompt flux. To minimize 
the statistical uncertainty the prompt flux has been averaged over the flavor (u^, v e and 
fi), since the cross-sections and branching ratios in the considered channels are nearly equal 
[62] . The results are shown in Fig. [8] In general, the absolute flux generated with QGSJet 
is at least factor of ~ 2 lower in the range compared to the two reference calculations, the 
spectral index shows compatible characteristics up to energies of hundreds of TeV. Since 
the crossover between the conventional and the prompt flux occurs in a region above the 
knee, where the spectral index of the conventional flux has its maximum, the uncertainty is 
too large for a detailed discussion. Therefore, we restrict ourselves to the discussion of the 
influence of the primary model. 

In the region below 100 TeV where the conventional fluxes dominate, the variation of the 
prompt component due to the primary model is smaller than the differences between the 
theoretical predictions in Fig. [8j Above the knee the influence of the cosmic ray intensity and 
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FIG. 8. (Prompt) Flux of leptons originating from charm particles, calculated with QGSJet-Olc. 

composition is evident. Primary models, which do not include a second Galactic or a third 
extragalactic component, suffer from a steep cutoff (poly-gonato and ZS). The differences 
between the cHGp and cHGm models show, that for the highest energies the modeling of an 
extragalactic component is crucial and that different composition scenarios of the primary 
flux influence the total rate of prompt particles. 



E. Estimation of the theoretical uncertainty 

The total theoretical uncertainty of the conventional flux of atmospheric muons and 



neutrinos, derived from the results of this calculation, is shown in Fig. [9] and Tab. VE The 
shaded bands were calculated using the average particle spectrum: 

(%(E)) = — 1— £ £ %(E, M, $ c ) (5) 

The upper and lower boundaries of the uncertainty bands are then given by 

O^IE) = max fL — , r . o (E) = mm fL — , . . . 6 

This is for the combined uncertainty due to the interaction model and the primary flux 
model. Since the GH and the poly-gonato model do not consider a steepening and the 
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FIG. 9. Total theoretical uncertainty derived from the variation of the interaction model is repre- 
sented by the inner solid bands. The shaded bands are derived from the variation of the interaction 
model using all available primary models (cHGp, cHGm, ZS/PAMELA for the medium shaded 
bands and cHGp, cHGm, ZS/PAMELA, GH and poly-gonato for the outer shaded bands). 



disagreement of single power law fits to observations of direct measurements, we provide in 
the medium shaded band the uncertainty when using cHGp, cHGm and ZS as candidate 
spectra only. This result suggests that up to 1 TeV neutrino or muon energy, the uncertainty 
of the flux is dominated by the hadronic interaction model, while the flux of primary cosmic 
rays and the composition (< 100 TeV/nucleus) are relatively well known. Above this energy 
the different assumptions about the origin and shape of the knee become dominant and 
result in additional 15% to the uncertainties of the hadronic interaction models. At lepton 
energies approaching 100 TeV the uncertainties due to the primary flux are decreasing, since 
the primary nuclei responsible for these particles are from an energy range at the knee, where 
the primary models are compatible with each other. 

The pure interaction model uncertainty was calculated by fixing the primary model $c 



to cHGp. The result shows the features discussed in Sec. V A The weaker dependence 



on the representation of kaons in the air showers and a better overall agreement of the 
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TABLE I. Theoretical uncertainties of the conventional atmospheric lepton fluxes, given in %. ALL 

- variation of the combinations between the all primary models (cHGp, cHGm, ZS-PAMELA, poly- 
gonato and GH) and all hadronic interaction models HIM (SIBYLL-2.1, QGSJetOlc, QGSJet-II). 
RECENT contains the same set as ALL, excluding the GH and the poly-gonato models. HIM only 

- primary model is fixed to cHGp and all interaction models are varied. The average values (av) 
are calculated using all data points from the bands. 



Type 


ALL + HIM 


RECENT + HIM 


HIM only 


Energy (TeV) 


0.1 


1 


10 


100 


av 


0.1 


1 


10 


100 


av 


0.1 


1 


10 


100 


av 




+32 
-23 


+42 
-27 


+47 
-30 


+39 
-41 


+42 
-32 


+29 
-17 


+32 
-19 


+35 
-27 


+29 
-23 


+32 
-22 


+26 
-16 


+30 
-16 


+26 
-18 


+27 
-20 


+27 
-20 




+23 
-18 


+35 
-25 


+39 
-26 


+33 
-39 


+35 
-30 


+19 
-13 


+25 
-16 


+28 
-23 


+24 
-16 


+25 
-19 


+17 
-12 


+23 
-14 


+19 
-13 


+18 
-13 


+20 
-13 




+12 
-15 


+20 
-20 


+29 
-22 


+23 
-30 


+24 
-24 


+10 
-13 


+12 
-11 


+19 
-19 


+13 
-10 


+15 
-13 


+7 
-9 


+11 

-9 


+11 

-8 


+11 
-7 


+11 

-7 



pion performance between the interaction models leads to the determination of the muon 
flux with a high precision (< 11 %) up to hundreds of TeV. The important role of kaons 
for atmospheric neutrino production results in high uncertainties for both neutrino flavors. 
For the production of muon neutrinos, pions and the i\"/7r-ratio play a bigger role, thus 
resulting in higher uncertainties for muon- than for electron neutrinos. Our results for the 
interaction model uncertainty are somewhat higher when compared to the detailed study 
of uncertainties in atmospheric neutrino fluxes by Barr et al. [87], which predicts for muon 
neutrinos at 1 TeV a total (gaussian) uncertainty of 30% due to hadronic interactions and 
of 40% if the primary flux is taken into account. 



1. Uncertainties due to the composition in the UHE regime 

Since the flux of muons and neutrinos at energies above hundreds of TeV is dominated 
by prompt muons and neutrinos, only QGSJet-Olc with an enabled charm component is 



suitable for studies. In Fig. 10, the same approach as in the previous section has been 



chosen to identify the uncertainty of the UHE component due to the uncertain composition 

26 



10 5 10 6 10 7 E / GeV 

FIG. 10. Theoretical uncertainty of the prompt atmospheric leptons due to the composition of the 
extragalactic component component. 

of extra-galactic cosmic rays, resulting in less than ±2% at 100 TeV and ±40 % at hundreds 
of PeV. 

VI. SUMMARY AND DISCUSSION 

We have developed a full Monte-Carlo calculation scheme which is capable to calculate 
muon neutrino, electron neutrino and muon fluxes up to 100 TeV, with a statical accuracy 
around a few percent. For surface energies up to 100 PeV, it was possible to assess the 
influence of the spectrum and composition of primary nuclei on the calculation. With 
respect to the increasing sensitivity of modern neutrino telescopes up to these energies, the 
influence of the knee of cosmic rays has been studied from the perspective of lepton spectra 
at the surface. 

The asymmetric uncertainties in the calculation of the fluxes have been assessed by car- 
rying out the calculation using several models of the primary cosmic ray flux and three 
interaction models, which individually represent different assumptions about hadronic inter- 
actions. It has been found that the uncertainties in the calculation of the atmospheric muon 
flux are significantly smaller, compared to the uncertainties of the atmospheric neutrino 
flux. As it has been shown in the study of flavor and charge ratios, this behavior can be 
explained by the insufficiently known contribution of kaons in the atmospheric cascade and 
the higher importance of kaons for the neutrino flux. 

In particular, uncertainties of the conventional atmospheric neutrino flux are important 
as a dominant source of systematic uncertainties in the search for astrophysical high-energy 
neutrinos. Recent searches for different astrophysical neutrino signals with IceCube are 
based on the Gaisser-Honda parametrization of the spectrum with an assumed systematic 
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uncertainty of 25% [7H [88] . Here, it could be shown that depending on the energy, the 
uncertainty of the atmospheric neutrino flux is ±22% 011 average when using a realistic 
cosmic ray spectrum. While this number is dominated by the interaction model, there is 
still some significant contribution from the primary flux models. If a better understanding 
of the composition and spectral behavior up to the knee and above can be achieved, the 
total uncertainty can be reduced further. In addition, it is expected that the inclusion of 
new LHC data can reduce the systematic uncertainties coming from the interaction models 
[81?] at the energy of the knee. 

Using the charm option provided with QGS Jet-Olc, it was possible to calculate the prompt 
component of atmospheric leptons. Although, the absolute value of the prompt flux is 
significantly lower than expected from other calculations, we were able to show the role of 
the primary flux model in this type of calculations. The prompt flux has been identified to 
be sensitive to the composition of the extragalactic cosmic rays. 
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